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The nonlinear gyrokinetic code GS2 has been extended to treat non-axisymmetric stellarator geometry. Elec- 
tromagnetic perturbations and multiple trapped particle regions are allowed. Here, linear, collisionless, elec- 
trostatic simulations of the quasi- axisymmetric, three-field period National Compact Stellarator Experiment 
(NCSX) design QAS3-C82 have been successfully benchmarked against the eigenvalue code FULL. Quanti- 
tatively, the linear stability calculations of GS2 and FULL agree to within ~ 10%. 



I. INTRODUCTION 

One of the most important issues for magnetic fu- 
sion is the confinement of heat and particles. Turbulent 
transport (most likely the result of drift wave instabil- 
ities) causes a significant amount of heat loss in toka- 
maks and spherical toriii Neoclassical transport, on the 
other hand, can often account for the poor confinement 
in traditional stellarators. 2 However, modern stellarator 
designs, such as Wendelstein 7-AS (W7-AS)j£ Wendel- 
stein 7-X (W7-X)^ the National Compact Stellarator 
Experiment (NCSX), 6 the Large Helical Device (LHD) 7 , 
and the Helically Symmetric Experiment (HSX)£~— have 
shown or are designed to have improved neoclassical con- 
finement and stability properties. Understanding plasma 
turbulence and transport could further improve the per- 
formance of stellarators. Progress in design of stellara- 
tors for optimal transport has been made by coupling 
the gyrokinetic code GENE 33 with the configuration op- 
timization code STELLOPT.— li 2 - 

Gyrokinetic studies of drift-wave-driven turbulence in 
stellarator geometry are relatively recent and comprehen- 
sive scans are scarce. Most of these studies were done 
using upgraded versions of well-established axisymmet- 
ric codes which include comprehensive kinetic dynam- 
ics (multispecies, collisions, finite beta) to the more gen- 
eral case of non-axisymmetric stellarator geometry, in the 
flux tube limit. The first non-axisymmetric linear gy- 
rokinetic stability studies, for both the ion-temperature- 
gradient-driven (ITG) mode and the trapped-electron 
mode (TEM), were done with the linear eigenvalue FULL 
codeji£~— including a comparison of stability in nine 
stellarator configurations, 16 Extensive studies have been 
done with the upgraded GENE code, including the first 
nonlinear gyrokinetic simulations. — More recently, the 
GKV-X code, which uses the adiabatic electron approx- 
imation, has been used to analyze linear ITG modes 
and zonal flows in LHD and nonlinear studies are in 
progress J£ 

For this purpose, the axisymmetric nonlinear microin- 
stability code GS2 19 has been extended to treat the more 
general case of non-axisymmetric stellarator geometry. 



GS2 contains a full (except that the equilibrium distri- 
bution function is taken to be a Maxwellian) implementa- 
tion of the 5-D Frieman and Chen nonlinear gyrokinetic 
equation in the flux tube limit j 1 ^^ with an efficient par- 
allelization for modern supercomputers^ It treats elec- 
trons and an arbitrary number of ion species on an equal 
footing, and includes trapped particles, electromagnetic 
perturbations, and a momentum-conserving pitch-angle- 
scattering collision operator. The extension of the code 
to non-axisymmetric geometry not only retains all of the 
above dynamics of the axisymmetric version, but also 
allows, most importantly, multiple trapped particle re- 
gions and multiple totally-trapped pitch angles at a given 
theta grid point. (By "totally-trapped," we mean parti- 
cles with such a small parallel velocity that they are lim- 
ited to one grid point at the bottom of a well.) Tokamaks 
only have one trapped particle region, but as stellarators 
can have many deep, narrow magnetic wells which can 
trap particles (though NCSX has only a single deep well, 
with other shallow wells, and is a bridge in configuration 
space between tokamaks and other stellarators). In or- 
der to treat the trapped particles accurately, one needs to 
resolve these wells sufficiently with high grid resolution. 
With the GS2 modifications, we allow for more flexible, 
decoupled pitch angle and parallel spatial grids, relative 
to the original GS2 algorithm which required every grid 
point (6j) along the field line to correspond exactly to 
the turning point of trapped pitch angle (A, = fi/E) grid 
points 
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Beyond these extensions, a GS2 stellarator simula- 
tion requires different geometry codes to build its in- 
put grids than standard tokamak runs. For these non- 
axisymmetric simulations, the geometrical coefficients 
are based on a VMEC S' -' 3D MHD equilibrium, which 
is transformed into Boozer coordinates^ by the TERP- 
SICHORE code. 24 From this equilibrium, the VVBAL 
code 2 ^ constructs data along a chosen field line nec- 
essary for the microinstability calculations: B = |B|, 
the VB drift, the curvature drift, and the metric coeffi- 
cients. While these extensions were used to study HSX 
plasmas^ here we verify the non-axisymmetric exten- 
sion of GS2 through comparisons with FULL on NCSX 
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plasmas. Good agreement between the GS2 code and 
the FULL code in the axisymmetric limit has been ex- 
tensively demonstrated previously . 27 ' 28 While the non- 
axisymmetric upgrade of GS2 retains the nonlinear dy- 
namics, in these studies we focus on systematic scans of 
gyrokinetic linear stability. 

The organization of this paper is as follows. The 
NCSX equilibrium used for the benchmark is described 
in Section II. Comparisons between the GS2 code and the 
FULL code in non-axisymmetric geometry over a range 
of parameters including r\ = L n / Lt (where L n is the 
density gradient scale length and Lt is the temperature 
gradient scale length), k y pi, Ti/T e , and geometrical coor- 
dinates are presented in Section III. Further results using 
the GS2 code to investigate effects of density and temper- 
ature gradients are presented in Section IV. Conclusions 
and a discussion of future work are given in Section V. 
Finally, Appendix A contains definitions of the normal- 
izations and radial coordinate used by GS2. 



II. THE QAS3-C82 EQUILIBRIUM 

All of the benchmark calculations use a VMEC equi- 
librium based on a 1999 NCSX design known as QAS3- 
C82^ which is shown in Figure 1. This configuration is 
quasi-axisymmetric with three field periods. It has an as- 
pect ratio of 3.5 and a major radius of 1.4 m. NCSX was 
designed to have good neoclassical transport and MHD 
stability properties and good drift trajectories similar to 
those in tokamaks. Strong axisymmetric components of 
shaping provide good ballooning stability properties at 
lower aspect ratio. Furthermore, the QAS3-C82 equi- 
librium has a monotonically increasing rotational trans- 
form profile which provides stability to neoclassical tear- 
ing modes across the entire cross section i 12 ' 29 

For most of these runs, we chose the surface at s = 
0.875 (s ~ ((r/a) 2 ) is the normalized toroidal flux) and 
the field line at a — ir/3 (a = £ — q9; ( is the Boozer 
toroidal angle, 9 is the Boozer poloidal angle). The cross- 
section at this point is the crescent shape, seen in Figure 
17 of Rcf. |30j. The coordinate along the field line is 8, 
the poloidal angle. At this surface, the safety factor q = 
2.118 and the average (3 (the ratio of the plasma pressure 
to the magnetic pressure) is {(i) = 0.01%. Lastly, the 
ballooning parameter 25 is 9q — 0, except in Figure 6. 

Figure 2 shows the variation of the magnitude of the 
magnetic field along a chosen magnetic field line. Reso- 
lution studies for the spatial grid used in the GS2 runs 
indicate that 330 theta grid points per poloidal period 
and about 90 pitch angles (A = fi/E) showed conver- 
gence in the growth rate to within 2%, however < 10% 
error is possible with coarser grids. It was also found 
that a 9 range extending from — 3ir to 37r was sufficient 
for a typical simulation grid, meaning that the eigenfunc- 
tions for the modes decayed to insignificant values before 
reaching these boundaries. (The endpoints of B(9) were 
increased slightly, by less than 1%, to be global maxima, 




FIG. 1. Equilibrium of NCSX design QAS3-C82 which is 
quasi-axisymmetric and has 3 field periods. 




(rad) 

FIG. 2. Standard B vs. 6 grid for QAS3-C82, with s = 0.875, 
a — 7r/3, and 6o = 0. 



per normal GS2 operations.) 

The equilibrium's geometry suggests unstable drift 
waves exist. The variations of (k±/n) 2 , where n is the 
toroidal mode number, and the curvature drift along the 
same chosen field line can be seen in Figures 3 and 4. By 
convention, positive curvature drifts are "bad" or desta- 
bilizing, while negative curvature drifts are "good" or 
stabilizing. Significant unstable modes occur where k± 
is small, which is near 9 = for this equilibrium, since 
instabilities are generally suppressed at large k± by FLR 
averaging. Also, because Figure 4 indicates that the cur- 
vature is bad in this region near 9 = 0, it is expected that 
unstable modes will appear here. 
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FIG. 3. Variation of (^) 2 (0) for QAS3-C82, with s = 0.875, 
a = 7r/3, and #o = 0. 




FIG. 4. Variation of the curvature drift frequency (w M = 
(k±/n) ■ b x [b ■ Vb]) (for n = 1) along 6> for QAS3-C82, with 
s = 0.875, a = tt/3, and 6 = 0. 



III. BENCHMARKS WITH FULL 

Comparisons between the GS2 code and the FULL 
code in non-axisymmetric geometry over a range of pa- 
rameters using the QAS3-C82 equilibrium show linear 
agreement for our standard case, whose local parameters 
are shown in Table I. The product of the perpendicu- 
lar wave number and the gyroradius at 8 — 0, k y pi, is 
0.3983 (where the toroidal mode number n = 25; see 
App. A) for all cases unless otherwise specified. The 
standard case is relatively close to the edge, which ac- 
counts for the low values of ion temperature, Tj, electron 
temperature, T e , and relatively large values for the gra- 
dients. The parameter, r\ = L„/L^ is usually 77 = 3, 
placing most of our studies in an ITG regime (see Figure 
7). Correspondingly, ctjv/L„i — / L ne — 13.096 and 
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TABLE I. The set of local parameters used in a standard case 
microinstability simulation based on the QAS3-C82 equilib- 
rium. Note: a at is not the minor radius; it is discussed in 
App. A. 



<in / Lti = cln I Ltc = 39.288. The major radius is ap- 
proximately R I Am. The normalizing scale length is 
ajy = n/k±_{6 = 0) = 0.352m, not the minor radius, and 
is described in detail in App. A. These studies are done 
with electrons and deuterium ions. 

Previously, FULL scans showed that the largest lin- 
ear growth rate occurs at flux surface label s = 0.875 
(corresponding to a minor radius of r/a ss -y/s ~ 0.94), 
for a = 7r/3 and #0 = 0. GS2 and FULL scans 
over a and #0 (Figures 5 and 6) adopted this s value. 
The toroidal mode number, n, was fixed at 25 (thus, 
kyPi = -j^-pi varied for each data point, because from 
App. A, = l/|Va| and pi oc 1/B a vary). These fig- 
ures indicate good agreement between the GS2 code and 
the FULL code. The maximum growth rate in Figure 
5 occurs for a = 7r/3, and GS2 and FULL agree well 
around this value. In Figure 6, GS2 and FULL again 
agree well around the growth rate peak at 8q = 0. 

In all further calculations presented in this paper, s — 
0.875, a = 7r/3 and 9q — 0, the location of the maximum 
growth rate. 

We used GS2 to find the instability growth rate depen- 
dence on 77 = L n /Lj> and compared it with FULL. The 
total pressure gradient was kept fixed to maintain con- 
sistency with the MHD equilibrium. Both codes found 
large growth rates at low 77 (high density gradient) and 
high 77 (high temperature gradient) (Figure 7), and agree 
well, though it can be seen in the frequencies that GS2 
found a mode switch earlier than FULL. This can happen 
since GS2 automatically finds the most unstable mode, 
whereas FULL usually finds the mode closest to the ini- 
tial guess provided to the root finder. In fact, there 
are three distinct eigenmodes within these regimes of 
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FIG. 6. (color online) Variation of 7 and u) r with So at 
constant s = 0.875 and a — n/3 with 77* = r/ e = 3 and 
fc H pi = 0.3983. 



FIG. 5. (color online) Variation of 7 and ui r with a at constant 
s — 0.875 and So = with rji — rj e = 3 and k y pi{a = = 
0.3983. 



77: at small 77, even-symmetry TEM modes dominate; at 
medium 77, odd-symmetry TEM modes dominate; and at 
larger values of 77, an even-symmetry ITG-driven mode 
dominates 15 (Figure 8). This is typical of an equivalent 
axisymmetric configuration. 31 

Benchmarks with FULL for scans over Ti/T e , shown in 
figure 9, were also successful. For this scan, T e was varied 
while Ti was kept constant at IkeV. As Ti/T e increases, 
at this very large value of R/L^ ~ 157, the linear growth 
rate falls slowly due, most likely, to an enhancement of 
shielding by adiabatic electrons at large ^/T,-/T e . This is 
a very well-known phenomenon in tokamaks. 

Comparison scans over k y pi for 77 = and 77 = 3 are 
shown in figure 10. For the 77 = curve, the dominating 
eigenmodes are even in the ranges 0.1 < k y pi < 0.2 and 
0.6 < kyPi < 1.1. Overall, the results from the GS2 code 
and the FULL code agree well; growth rates differ by at 
most ~ 10% except at transitions between modes. 

We found high frequency, electron-temperature- 
gradient-driven (ETG) modes with GS2 at short wave- 
lengths (Figure 11) in the extended k y pi spectrum for the 
case of f] = 3. This was not checked with FULL. 
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FIG. 7. (color online) Variation of 7 and uj r with 77^ = ?7 e with 
k y pi = 0.3983. 
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FIG. 8. (color online) Variation of the normalized GS2 eigen- 
functions of electrostatic, collisionless toroidal drift modes 
along the field line at rj = 3 (top figure) and at rj = 0.5 
(bottom figure) with k y pi = 0.3983. 
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FIG. 10. (color online) Variation of 7 and u) r with k y pi. Cir- 
cles: GS2, 77 = 0; triangles: FULL, rj = 0; squares: GS2, 
rj = 3; diamonds: FULL, 77 = 3. 
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FIG. 9. (color online) Variation of 7 and uj r with Ti/T e with 
k y pi — 0.3983 and rji = rj e = 3. 



FIG. 11. (color online) Extended variation from GS2 of 7 and 
u> r with fe^pi and rji = r\ & = 3. 
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FIG. 13. (color online) Variation from GS2 of 7 and u) r 
with aN/Ln with fc^pi = 0.3983. Circles: ajv/L„ = 
52.4, a N /L Te = 0.0; triangles: a N /L n = 13.1, a N /L Te = 
39.3; squares: ajv/L n = 0.0, ajv/Z/Te = 0.0. 



FIG. 12. (color online) Variation of 7 and uj t with ajv /L n with 



= 0.3983. Circles: on/Lt 
39.3; squares: ajv/Lr = 44.9. 
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IV. CRITICAL GRADIENTS FOR LINEAR 
INSTABILITY 

GS2 was also used to search for critical density gradi- 
ents and temperature gradients; i.e. to see whether gra- 
dients exist at which all drift wave modes are stabilized. 
Note that for the next series of figures, the normalizing 
length for the density and temperature gradient length 
scales is defined as = (n/k±)(6 = 0) ~ 0.352m (see 
App. A). 

Figure 12 shows a scan over the density gradient at 
various ion and electron temperature gradients. The re- 
sults are inconsistent with the equilibrium pressure gra- 
dient, as the density gradient was increased at constant 
temperature gradient. However, because the equilibrium 
beta is so small (~ 0.01%), the effect of the variation 
of the pressure gradient is negligible. We see that there 
is no nonzero critical density gradient threshold, even in 
the absence of temperature gradients. There are switches 
in eigenmode symmetry from even to odd as <in / L n in- 
creases, or all cin/Lt values. 

However, a critical ion temperature gradient for an 
ITG-driven mode was found at a^/Lri ~ 2 (or = 
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FIG. 14. (color online) Variation from GS2 of 7 and u r with 
aN/LTe for the case of Fig. 2 with k y pi — 0.3983. Cir- 
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_R_ on. ~ /±±n_ — i n the absence of all other gradients 

(Figure 13). Likewise, a critical electron temperature 
gradient for a TEM-driven mode was found at j^- ~ 2 

in the absence of all other gradients (Figure 14). 

V. CONCLUSIONS 

The nonlinear gyrokinetic code GS2 has been extended 
to treat non-axisymmetric stellarator geometry. Geomet- 
ric quantities required for the gyrokinetic simulations are 
calculated from a VMEC-generated equilibrium using the 
VVBAL code and are further described in App. A. 

Linear, collisionless, electrostatic simulations of the 
quasi- axisymmetric, three-field period NCSX stellarator 
design QAS3-C82 have been successfully benchmarked 
with the eigenvalue code FULL for scans over a range of 
parameters including 77, k y pi, Ti/T e , a, and 9q. Quan- 
titatively, the linear stability calculations of GS2 and 
FULL agree to within about 10% of the mean, except 
at transitions between modes. Further results using only 
GS2 included short wavelength modes, odd parity, faster 
growing modes, and the effect of individual density and 
temperature gradients. 

Future work will include the exploration of the effects 
of collisionality and electromagnetic dynamics, investiga- 
tion of finite beta equilibria, and, most significantly, the 
effects of nonlinear dynamics. A benchmark of stellarator 
studies is underway between GS2 and the continuum gy- 
rokinetic code GENE22 for NCSX, as well as stellarators 
W7-AS and W7-X. 

GIST— is now capable of creating GS2 geometry data 
files and will be used in the future, along with the new 
GS2 grid generator. 
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Appendix A: Geometry Details 

In order to make the simulation grid for these GS2 stel- 
larator runs, VMEC creates the 3-D MHD equilibrium, 
TERPSICHORE transforms it into Boozer coordinates, 
and VVBAL calculates necessary geometric coefficients 



along a specified field line. Then, GS2's grid generator, 
Rungridgen, creates the final grid for use in the microin- 
stability calculations. (A new grid generator is in pro- 
duction, which will be used for further GS2 stellarator 
calculations.) The normalizations of geometric quantities 
change between these codes, and knowing them in detail 
is required for benchmarks between gyrokinetic codes. 
We define the normalizing length, a^, in App. A. 2. 

In GS2, the field-aligned coordinate system is (p, a, 9). 
9 is the poloidal angle and distance along the field line. 
The magnetic field takes the form B = Vq x VP, where 
a = £ — q9 is the field line label. The radial coordinate, 
p, can differ between codes, and we define it in App. A.l. 
More details of general geometry for GS2 are documented 
in App. A of Ref. [M 

1. Radial coordinate, p 

VMEC and TERPSICHORE use the normalized 
toroidal flux surface label s = $/& e dge ~ (( r / a ) 2 ) as 
the radial coordinate, p. In the customized version of 
VVBAL used here, the radial coordinate is transformed 
to p = "i>N = ^ / i a %B a ), where "J at is the normalized 
poloidal flux. 

Because Rungridgen uses VVBAL output without 
modification, here dp/dipN = 1- (In Ref. l35l . the def- 
inition of the geometry coefficients include the variable 
dp/dipN, which can be used to choose the radial coordi- 
nate.) 

2. Normalizing Quantities, B a and ajv 

The normalizing magnetic field is B a = (B), where 
(B) is a theta-average, not weighted to be a flux-surface 
average (Ref. [HI chooses B a differently) . 

The normalizing length is a^, given for these calcula- 
tions by VVBAL as 

a N = , = 7^— r- (Al) 

V|kj 2 (0 = O,0o = O) |Va| 

GS2 treats perturbed quantities as A = A(9)exp(iS), 
where k± = VS = nV(a + q9 ) = riV[C - q(9 - 6 )]; 
n is the toroidal mode number. (In non-axisymmetric 
devices, n is not a conserved quantum number, because 
toroidal variations in the equilibrium give coupling be- 
tween n modes. However, in the small-p*, high-rt limit, 
this coupling is weak, and n can just be considered a 
coefficient to select a particular value of k±.) 

In the notation of Eqn. A. 11 of App. A in Ref. |H, 

|k ± | 2 = \VS\ 2 = k 2 9 \g x + 29 Q g 2 + 9 2 Q g 3 \ (A2) 

where g±, and 53 are coefficients in the geometry 
file written by VVBAL and read by GS2. Also, kg = 
k y = n/aM- (The GS2 variable aky is defined 
with p t oc 1/B a .) 



7 



In the notation of Eqn. 7 of Ref. 



|k,| 2 = 



*' 2 (s) 



[C p + C s (9-e )+C q (9-9o) 2 }, (A3) 



where ^/g is the Jacobian, C p , C s , and C q are defined in 
section II of Ref. 
So, VVBAL writes: 

.91 = a 2 N [C p + C s e + C q 9 2 ] (A4) 
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